﻿//Class-filter of signal
namespace Diplom.Wavelets
{
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using Diplom.NumericArrays;

    
    

    public static class Filter
    {
        #region Types
        private struct Decomposition
        {
            public int WaveletOrder;
            public int SignalLength;
            public Numeric1DArray Approx;
            public Numeric2DArray Details;
        }

        private enum FiltrationResult
        {
            Success = 0,
            WrongOrder = 1,
            SignalTooShort = 2,
            WrongLength = 3
        }
        #endregion

        #region Constants
        private const int NWavelets = 10;

        private static double[][] Daubeshies = new double[][]{
                                                 new double[]{1},
                                                 new double[]{0.5,0.5},
                                                 new double[]{0.34150635094622,0.59150635094587, 0.15849364905378,-0.09150635094587},
                                                 new double[]{0.23523360389270,  0.57055845791731,  0.32518250026371, -0.09546720778426,
                                                    -0.06041610415535,  0.02490874986589 },
                                                 new double[]{0.16290171402562,  0.50547285754565,  0.44610006912319, -0.01978751311791,
                                                    -0.13225358368437,  0.02180815023739,  0.02325180053556, -0.00749349466513},
                                                 new double[]{0.11320949129173,  0.42697177135271,  0.51216347213016,  0.09788348067375,
                                                    -0.17132835769133, -0.02280056594205,  0.05485132932108, -0.00441340005433,
                                                    -0.00889593505093,  0.00235871396920},
                                                 new double[]{0.07887121600143,  0.34975190703757,  0.53113187994121,  0.22291566146505,
                                                    -0.15999329944587, -0.09175903203003,  0.06894404648720,  0.01946160485396,
                                                    -0.02233187416548,  0.00039162557603,  0.00337803118151, -0.00076176690258},
                                                 new double[]{0.05504971537285,  0.28039564181304,  0.51557424581833,  0.33218624110566,
                                                    -0.10175691123173, -0.15841750564054,  0.05042323250485,  0.05700172257986,
                                                    -0.02689122629486, -0.01171997078235,  0.00887489618962,  0.00030375749776,
                                                    -0.00127395235906,  0.00025011342658},
                                                 new double[]{0.03847781105406,  0.22123362357624,  0.47774307521438,  0.41390826621166,
                                                    -0.01119286766665, -0.20082931639111,  0.00033409704628,  0.09103817842345,
                                                    -0.01228195052300, -0.03117510332533,  0.00988607964808,  0.00618442240954,
                                                    -0.00344385962813, -0.00027700227421,  0.00047761485533, -0.00008306863060 },
                                                 new double[]{0.02692517479416,  0.17241715192471,  0.42767453217028,  0.46477285717278,
                                                    0.09418477475112, -0.20737588089628, -0.06847677451090,  0.10503417113714,
                                                    0.02172633772990, -0.04782363205882,  0.00017744640673,  0.01581208292614,
                                                    -0.00333981011324, -0.00302748028715,  0.00130648364018,  0.00016290733601,
                                                    -0.00017816487955,  0.00002782275679 },
                                                 new double[]{0.01885857879640,  0.13306109139687,  0.37278753574266,  0.48681405536610,
                                                    0.19881887088440, -0.17666810089647, -0.13855493935993,  0.09006372426666,
                                                    0.06580149355070, -0.05048328559801, -0.02082962404385,  0.02348490704841,
                                                    0.00255021848393, -0.00758950116768,  0.00098666268244,  0.00140884329496,
                                                    -0.00048497391996, -0.00008235450295,  0.00006617718320, -0.00000937920789}
                                             };
        #endregion

        #region metods

        public static void Denoise(ref Numeric1DArray array,int order,int level)
        {
            //Init();
            Decomposition decomp;
            if (FiltrationResult.Success!=Decompose(ref array, order, level,out decomp)) return;
            for (int i = 0; i < decomp.Details.NSize; i++)
            {
                for (int j = 0; j < decomp.Details.Items[i].Length; j++)
                {
                    decomp.Details.Items[i][j] = 0;
                }
            }
            Reconstruct(ref decomp, out array);
            array.InitCharacters();
        }
        
        /// <summary>
        /// Signal Decomposition
        /// </summary>
        /// <param name="signal">Signal</param>
        /// <param name="Order">Wavelet Order</param>
        /// <param name="Level">Level of decomposition</param>
        /// <param name="decomposition">Signal Decomposition</param>
        /// <returns>Result</returns>
        private static FiltrationResult Decompose(ref Numeric1DArray signal, int Order, int Level, out Decomposition decomposition)
        {
            Decomposition decomp = new Decomposition();
            if (Order < 1 || Order > NWavelets)
            {
                decomposition = decomp;
                return FiltrationResult.WrongOrder;
            }
            Numeric1DArray Filter = new Numeric1DArray(Daubeshies[Order].Length);
            for (int i = 0; i < Filter.Length; i++)
            {
                Filter.Insert(i, Daubeshies[Order][i]*Math.Sqrt(2.0));
            }

            decomp.WaveletOrder = Filter.Length / 2;
            decomp.SignalLength = signal.Length;
            decomp.Approx = new Numeric1DArray(signal.Length);
            for (int i = 0; i < signal.Length; i++)
            {
                decomp.Approx.Insert(i, signal.Items[i]);
            }
            decomp.Details = new Numeric2DArray(Level);
            for (int i = 0; i < Level; i++)
            {
                FiltrationResult result = DWT(ref decomp.Approx, ref Filter, ref  decomp.Details.Items[i]);
                if (result != FiltrationResult.Success)
                {
                    decomposition = decomp;
                    return result;
                }
            }
            
            decomposition = decomp;
            return FiltrationResult.Success;
        }

        /// <summary>
        /// Forward Discrete Wavelet Transform
        /// </summary>
        /// <param name="Signal">Signal</param>
        /// <param name="Filter">Filter</param>
        /// <param name="Details">Details of Transform</param>
        /// <returns>Result</returns>
        private static FiltrationResult DWT(ref Numeric1DArray Signal, ref Numeric1DArray Filter,ref Double[] Details)
        {
            Numeric1DArray Buffer;
            int K = Signal.Length;
            int L = Filter.Length;
            int N = K;
            if ((K & 1) == 1) N++;
            if (L > N) return FiltrationResult.SignalTooShort;
            Buffer = new Numeric1DArray(N + L - 2);
            for (int i = 0; i < K; i++)
            {
                Buffer.Insert(i, Signal.Items[i]);
            }
            if ((K & 1) == 1) Buffer.Insert(K, (Signal.Items[0] + Signal.Items[K - 1]) / 2);
            if (L > 2)
            {
                for (int i = 0; i < L - 2; i++)
                {
                    Buffer.Insert(N + i, Signal.Items[i]);
                }
            }
            N /= 2;
            K = L / 2;

            Signal = new Numeric1DArray(N);
            Details = new double[N];
            for (int i = 0; i < N; i++)
            {
                int i0 = i * 2;
                Signal.Items[i] = 0;
                Details[i] = 0;
                for (int j = 0; j < K; j++)
                {
                    int j0 = j * 2;
                    int j1 = j0 + 1;
                    int j2 = j0 + 2;
                    Signal.Items[i] += (Buffer.Items[i0 + j0] * Filter.Items[j0] + Buffer.Items[i0 + j1] * Filter.Items[j1]);
                    Details[i] += ((Buffer.Items[i0 + j0] * Filter.Items[L - j1] + Buffer.Items[i0 + j1] * Filter.Items[L - j2]));
                }
            }
            return FiltrationResult.Success;
        }

        /// <summary>
        /// Signal reconstruction
        /// </summary>
        /// <param name="decomp">Decomposition</param>
        /// <param name="Signal">Signal</param>
        /// <returns>Result</returns>
        private static FiltrationResult Reconstruct(ref Decomposition decomp, out Numeric1DArray Signal)
        {
            Numeric1DArray signal = decomp.Approx;
            if (decomp.WaveletOrder < 1 || decomp.WaveletOrder > NWavelets)
            {
                Signal = signal;
                return FiltrationResult.WrongOrder;
            }
            Numeric1DArray Filter = new Numeric1DArray(Daubeshies[decomp.WaveletOrder].Length);
            for (int i = 0; i < Filter.Length; i++)
            {
                Filter.Insert(i, Daubeshies[decomp.WaveletOrder][i]*Math.Sqrt(2.0));
            }
            for (int i = decomp.Details.NSize-1; i >= 0; i--)
            {
                bool isOdd = (((i > 0) && ((decomp.Details.Items[i - 1].Length % 2) == 1)) || ((i == 0) && ((decomp.SignalLength % 2) == 1)));
                signal = IDWT(ref Filter, ref decomp.Approx, ref decomp.Details.Items[i], isOdd);
                decomp.Approx = new Numeric1DArray(signal.Length);
                for (int j = 0; j < signal.Length; j++)
                {
                    decomp.Approx.Insert(j, signal.Items[j]);
                }
            }
            Signal = signal;
            if (signal.Length != decomp.SignalLength)
            {
                return FiltrationResult.WrongLength;
            }

            return FiltrationResult.Success;
        }

        /// <summary>
        /// Reverse Discrete Wavelet Transform
        /// </summary>
        /// <param name="Filter">Filter</param>
        /// <param name="Approx">Approximation</param>
        /// <param name="Details">Detail</param>
        /// <param name="isOdd">isOdd flag</param>
        /// <returns>Signal</returns>
        private static Numeric1DArray IDWT(ref Numeric1DArray Filter, ref Numeric1DArray Approx, ref Double[] Details, bool isOdd)
        {
            Numeric1DArray result = null;
            Numeric1DArray AppBuf = null;
            Numeric1DArray DetBuf = null;
            int Na = Approx.Length;
            int Nd = Details.Length;
            int L = Filter.Length;
            int K = (L / 2) - 1;
            int N = Math.Max(Na, Nd);
            AppBuf = FillBuff(ref Approx, N, Na, K);
            DetBuf = FillBuff(ref Details, N, Nd, K);
            double[] res = new double[2 * N];
            for (int i = 0; i < N; i++)
            {
                int i0 = i * 2;
                int i1 = i0 + 1;
                res[i0] = res[i1] = 0;
                for (int j = 0; j <= K; j++)
                {
                    int j0 = j * 2;
                    int j1 = j0 + 1;
                    int j2 = j0 + 2;
                    res[i0] += (AppBuf.Items[i + j] * Filter.Items[L - j2] + DetBuf.Items[i + j] * Filter.Items[j1]);
                    res[i1] += (AppBuf.Items[i + j] * Filter.Items[L - j1] + DetBuf.Items[i + j] * Filter.Items[j0]);
                }
            }
            if (isOdd)
            {
                result = new Numeric1DArray(2 * N - 1);
            }
            else
            {
                result = new Numeric1DArray(2 * N);
            }
            for (int i = 0; i < result.Length; i++)
            {
                result.Insert(i, res[i]);
            }
            return result;
        }

        /// <summary>
        /// FillBuffer
        /// </summary>
        /// <param name="source">Source Array</param>        
        /// <returns>Resut array</returns>
        private static Numeric1DArray FillBuff(ref double[] source, int N, int Nb, int L)
        {
            Numeric1DArray Result = new Numeric1DArray(N + L);
            if (Nb > 0)
            {
                for (int i = 0; i < Nb; i++)
                {
                    Result.Insert(L + i, source[i]);
                }
            }
            if (Nb < N)
            {
                for (int i = 0; i < N - Nb; i++)
                {
                    Result.Items[L + Nb + i] = 0;
                }
            }
            if (L > 0)
            {
                for (int i = 0; i < L; i++)
                {
                    Result.Insert(i, Result.Items[N+i]);
                }
            }
            return Result;
        }

        /// <summary>
        /// FillBuffer
        /// </summary>
        /// <param name="source">Source Array</param>        
        /// <returns>Resut array</returns>
        private static Numeric1DArray FillBuff(ref Numeric1DArray source, int N, int Nb, int L)
        {
            Numeric1DArray Result = new Numeric1DArray(N + L);
            if (Nb > 0)
            {
                for (int i = 0; i < Nb; i++)
                {
                    Result.Insert(L + i, source.Items[i]);
                }
            }
            if (Nb < N)
            {
                for (int i = 0; i < N - Nb; i++)
                {
                    Result.Items[L + Nb + i] = 0;
                }
            }
            if (L > 0)
            {
                for (int i = 0; i < L; i++)
                {
                    Result.Insert(i, Result.Items[N + i]);
                }
            }
            return Result;
        }

        /// <summary>
        /// Simple denoise function
        /// </summary>
        /// <param name="source"></param>
        public static void SimpleDenoise(ref Numeric2DArray source, int length, double dispersion)
        {
            double sigma = Math.Sqrt(2 * Math.Log(length + 0.0));
            sigma *= dispersion;
            for (int i = 0; i < source.NSize; i++)
            {
                for (int j = 0; j < source.MSize; j++)
                {
                    if (Math.Abs(source.Items[i][j]) < sigma)
                    {
                        source.Insert(i, j, 0.0);
                    }
                    else
                    {
                        if (source.Items[i][j] > 0)
                        {
                            source.Insert(i, j, source.Items[i][j] - sigma);
                        }
                        else
                        {
                            source.Insert(i, j, source.Items[i][j] + sigma);
                        }
                    }
                }
            }
        }

        #endregion

    }
}
